function simulateROLs!(cc,lvnew,theta,temps)
    simulateRCs!(theta,cc,lvnew) #RC's
    simulateU0!(theta,cc,lvnew) #U0
    simulateU!(theta,cc,lvnew,temps) #U
end

function compareAlgorithms(cc,lvnew,cf_constants,cf_temps,theta,full=true)
    II = size(cf_temps.Ufull,2)
    Jtilde = size(cf_temps.Ufull,1)
    if full
        platform = trues(Jtilde-1)
        cap = [max.(cf_constants.seats, [sum(cc.enroll.==jj) for jj=1:size(lvnew.U,1)]); II]
    else
        platform = cc.platform
        cap = [cf_constants.seats+cf_constants.sobrecupo; 0]
    end
    # run CPDA
    setupInitialDA!(cf_temps,cf_constants,lvnew,platform,cc,theta)
    runDA!(cf_temps,cap)
    placement = copy(cf_temps.point_at)
    # run SPDA
    setupInitialDA!(cf_temps,cf_constants,lvnew,platform,cc,theta)
    runSPDA!(cf_constants,cf_temps,cap)
    placementSPDA = copy(cf_temps.point_at)
    return placement, placementSPDA
end


function simulateEverything!(dfs,cfctlname,theta,cc,lv,temps,cf_constants,cf_temps)
    II = size(lv.U,2)
    cap0 = [cf_constants.seats+cf_constants.sobrecupo; 0] #baseline, initial match
    cap0_after = [cf_constants.seats; II]
    lvnew = deepcopy(lv)
    simulateEverything!(dfs,cfctlname,theta,cc,lvnew,temps,cf_constants,cf_temps,cap0,cap0_after)
end

function simulateEverything!(dfs,cfctlname,theta,cc,lvnew,temps,cf_constants,cf_temps,cap0,cap0_after)
    simulateRCs!(theta,cc,lvnew) #RC's
    simulateU0!(theta,cc,lvnew) #U0
    simulateU!(theta,cc,lvnew,temps) #U
    # A, placement, enrollment, H, graduation computed inside "two-stage match"
    II = size(lvnew.U,2)
    s2013 = twoStageMatch(dfs,cfctlname,lvnew,cap0,cap0_after,theta.alpha,cc.platform,cf_constants,cf_temps,cc,theta, storeInitialPlacement=true)
end

function simulateRCs!(theta,ccnew,lvnew)
    (nrc,II) = size(lvnew.rc)
    for ii=1:II
        typ = ccnew.typ[ii]
        lvnew.rc[:,ii] .= rand(MvNormal(zeros(nrc),theta.Sigma_rc[:,:,typ]))
    end
end

function simulateU0!(theta,ccnew,lvnew)
    II = size(lvnew.U0,2)
    for ii=1:II
        typ = ccnew.typ[ii]
        lvnew.U0[1,ii] = view(ccnew.Xoo,ii,:)'*theta.betaOO0[:,typ] + rand(Normal(0.0,sqrt(theta.sigsqOO0[typ])))
        lvnew.U0[2,ii] = view(ccnew.Xoo,ii,:)'*theta.betaOO1[:,typ] + rand(Normal(0.0,sqrt(theta.sigsqOO1[typ])))
    end
end

function simulateU!(theta,ccnew,lvnew,temps)
    precomputeMeanUtilities!(theta,lvnew,ccnew,temps)
    II = size(lvnew.U0,2)
    for ii=1:II
        updatemu!(theta,lvnew,ccnew,temps,ii)
        lvnew.U[:,ii] .= temps.mu[:,ii] .+ rand.(Normal())
    end
end

#
# ccnew = deepcopy(dataset[2012][1])
# lvnew = dataset[2012][2]
# temps = dataset[2012][3]

# df_individuals = DataFrame(ind=1:II, type=cc.typ)
# df_means = DataFrame(inds=["All","Male Private","Male Public","Female Private","Female Public"])
# df_byIter = Dict{Symbol,Array{Float64,1}}()
# dfs = (df_byIter,df_means,df_individuals)
